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Examples  of  simple  zone  fire  models  are  analyzed.  These  models  illustrate  the 
nature  of  the  numerical  problems  commonly  encountered  in  zone  models  of  en- 
closure fires.  Often  these  difficulties  arise  in  the  solution  of  the  equations  for 
the  pressure  in  connected  rooms  because  the  pressure  equilibrates  much  more 
rapidly  than  other  dynamical  variables.  Since  these  models  are  very  simple, 
analytical  techniques  can  be  applied  and  some  insight  gained  regarding  the 
nature  of  these  problems.  The  models  consist  of  ordinary  differential  equa- 
tions coupled  with  algebraic  equations.  Singular  perturbation  methods  and 
phase  plane  analyses,  together  with  numerical  integration  of  the  appropriately 
nondimensionalized  equations,  are  employed  to  examine  the  stiff  nature  of  the 
equations  associated  with  these  models.  We  conclude  that  many  of  the  difficul- 
ties associated  with  numerical  integration  of  zone  fire  models  in  general  may 
be  circumvented  by  appropriate  analysis  of  the  zone  fire  model  equations. 


1 Introduction 

There  is  a long  history  of  analysis  of  the  dynamical  behavior  of  fires  in  buildings 
using  mathematical  models.  The  reason  for  development  of  the  mathematical  models 
and  their  use  in  practice  has  been  reviewed  in  [1].  The  original  mathematical  model 
of  a plume  used  in  zone  fire  models  was  developed  by  Morton,  Taylor  and  Turner 
[2].  Other  early  work  contributing  to  the  basic  development  of  these  models  includes 
experimental  [3]  and  theoretical  [4]  studies  of  the  effects  of  flow  through  openings 
induced  by  fires  in  enclosures;  analytical  examples  of  the  development  of  a stratified 
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ceiling  layer  and  the  filling  of  an  enclosure  by  the  heated  gases  [5];  analytical  exam- 
ples of  two  layer  modeling  of  the  smoke  movement  in  two-room  structures  [6];  and 
theoretical  study  of  the  flow  of  smoke  and  hot  gases  through  vents  [7]. 

In  this  brief  introduction,  we  do  not  attempt  a complete  review  of  the  literature 
on  zone  fire  models,  but  rather  try  to  give  a perspective  on  where  the  analysis  pre- 
sented contributes  to  these  efforts.  Mathematical  models  of  fires  have  commonly  been 
divided  into  two  categories,  field  models  and  zone  models.  A description  of  each  of 
these  types  of  model  as  well  as  their  relative  advantages  and  disadvantages  has  also 
been  discussed  in  [1].  A careful  derivation  of  the  equations  used  in  zone  fire  models 
from  the  basic  conservation  laws  utilizing  a control- volume  approach  has  been  given 
in  [8].  More  recently,  Forney  and  Moss  [9]  have  discussed  some  of  the  difficulties  long 
encountered  in  trying  to  integrate  numerically  the  zone  fire  models  for  general  prob- 
lems of  interest  in  the  fire  community.  Since  zone  models  have  been  so  useful  to  fire 
engineers  and  researchers,  and  since  numerical  integration  of  the  various  equations 
used  in  these  models  has  persistently  been  plagued  by  difficulties,  the  current  research 
is  presented. 

Examples  of  simple  zone  models  are  analyzed.  The  first  of  these  models  is  for  the 
pressure  in  a single  room  heated  by  a fire  and  vented  to  the  outside  by  a small  leak. 
This  model  has  been  solved  earlier  [5]  [9]  [10]  and  provides  the  basis  for  analyzing 
more  complex  models.  Here  the  single  ordinary  differential  equation  (ODE)  for  this 
model  is  made  dimensionless,  and  the  nondimensionalization  is  found  to  be  helpful 
for  analysis  of  more  complex  models. 

The  second  model,  the  coupled  pressure  equations  for  two  rooms  connected  to  each 
other  and  to  the  outside,  was  derived  in  the  report  of  Forney  and  Moss  [9]  because 
it  illustrates  in  very  simple  form  the  nature  of  the  problem  encountered  commonly 
in  zone  models  of  enclosure  fires.  Here,  the  equations  for  this  two-room  model  are 
made  dimensionless,  and  we  analyze  them  both  by  asymptotic  and  by  phase-plane 
methods.  The  methods  provide  insight  into  the  nature  of  the  numerical  problems 
commonly  encountered  in  zone  models. 

Two  appendices  are  also  included.  In  the  first,  we  provide  a derivation  of  a set 
of  differential  equations  which  govern  all  of  the  dependent  variables  (room  pressure, 
layer  height,  layer  temperatures  and  densities,  etc.)  for  a simple  two-layer  zone-fire 
model.  In  contrast  to  the  examples  given  in  the  main  body  of  this  report,  this  one 
demonstrates  that  the  equations  already  exhibit  stiffness  when  the  pressure  equation 
is  coupled  to  the  equations  governing  the  other  variables;  the  reason  for  this  stiffness 
is  that  the  pressure  generally  equilibrates  rapidly  compared  to  the  rate  at  which  the 
other  variables,  e.g.,  layer  height,  change  in  an  enclosure  fire.  In  the  second  appendix, 
we  present  the  commands  in  Mathematica  [11]  used  to  perform  the  calculations  and 
to  produce  the  figures  presented  in  this  report. 
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2 Conservation  Equations 

The  differential  equation  for  the  pressure  in  a simple  zone  fire  model  is  derived  in  this 
section  using  the  laws  of  conservation  of  mass  and  energy  together  with  the  equation 
of  state  for  an  ideal  gas.  Differential  equations  for  other  quantities  found  in  a zone  fire 
model  such  as  temperature  or  density  are  derived  in  the  general  case  in  [8]  [9];  a set 
of  differential  equations  for  a specific,  simple  two-layer  model  in  a single  enclosure  are 
derived  in  Appendix  A.  A zone  may  consist  of  a number  of  interior  regions  (usually 
an  upper  or  a lower  gas  layer).  The  basic  assumption  of  a zone  fire  model  is  that 
properties  such  as  temperatures  can  be  approximated  a^  uniform  throughout  the  zone. 
It  is  remarkable  that  this  assumption  seems  to  hold  for  as  few  as  two  gas  layers,  which 
is  the  model  considered  in  this  paper. 

Many  differential  equation  formulations  based  upon  these  assumptions  can  be  de- 
rived. One  formulation  can  be  converted  into  another  using  definitions  of  density, 
internal  energy  and  the  ideal  gas  law.  One  property  that  many  of  these  formulations 
share  is  the  presence  of  multiple  time  scales.  Physically,  the  pressure  in  a compart- 
ment equilibrates  much  quicker  than  densities  and  temperatures,  see  Appendix  A. 
Numerically,  this  property  is  known  as  stiffness  and  in  general  requires  the  use  of 
special  differential  equation  solvers  to  generate  efficient  solutions.  The  main  focus  of 
this  paper  then  is  to  show  how  an  appropriate  nondimensionalization  for  the  pressure 
equations  in  a zone  fire  model  can  be  used  with  analytic  techniques  such  as  phase- 
plane  analysis  and  singular  perturbation  methods  to  expose  and  exploit  the  presence 
of  these  multiple  time  scales. 

Each  differential  formulation  can  be  expressed  in  terms  of  mass  and  enthalpy  flow 
rates  denoted  rhu,  tni^  qu  and  qi  where  the  subscripts  L and  U refer  to  the  lower 
and  upper  layer  respectively.  These  flow  rates  represent  the  net  exchange  of  mass  or 
energy  between  zones  due  to  physical  phenomena  or  sub-models  such  as  fire  plumes, 
natural  and  forced  vents,  convective,  radiative  heat  transfer,  etc.  For  example,  a 
vent  exchanges  mass  and  energy  between  zones  in  connected  rooms,  a fire  plume 
typically  adds  heat  to  the  upper  layer  and  transfers  entrained  mass  and  energy  from 
the  lower  to  the  upper  layer,  and  convection  transfers  energy  from  the  gas  layers  to 
the  surrounding  walls. 

As  illustrated  in  Figure  1,  a compartment  can  be  divided  into  two  control  volumes, 
an  upper  layer  of  hot  gases  and  smoke,  and  a lower  layer  of  air.  The  fire  produces 
a plume  and  acts  as  a pump  to  transfer  mass  from  the  lower  to  the  upper  layer, 
adding  energy  to  the  transferred  fluid.  The  two  layer  model  is  quite  adequate  for 
many  applications  because  upper  and  lower  layers  as  described  are  often  observed 
experimentally  in  room  fires.  The  gas  in  each  layer  has  attributes  of  mass,  internal 
energy,  density,  temperature,  and  volume  denoted  respectively  by  m,,  F,,  p,,  T,,  and 
Vi  where  i = L for  the  lower  layer  and  i = U for  the  upper  layer.  The  compartment  as 
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a whole  has  the  attribute  of  pressure  P.  These  eleven  variables  are  related  by  means 
of  the  following  seven  constraints 


Tfl ' 

p,  = -^  (density)  = L,U  (1) 

Ei  = CuruiTi  (internal  energy)  ,i  = L,U  (2) 

P = RpiTj  (ideal  gas  law)  ,i  = L,U  (3) 

V = Vl  Vjj  (total  volume)  . (4) 


The  specific  heats  at  constant  volume  and  at  constant  pressure,  and  Cp,  the  uni- 
versal gas  constant,  R,  and  the  ratio  of  specific  heats,  7,  are  related  by 

7 = — , 

Cy 

R Cp  Cy  . 


The  first  law  of  thermodynamics  states  that  the  rate  of  increase  of  layer  internal 
energy  plus  the  rate  at  which  the  layer  does  work  by  expansion  is  equal  to  the  rate 
at  which  enthalpy  is  added  to  the  gas  (where  we  consider  the  enthalpy  added  as  that 
from  any  sources  minus  losses  to  the  walls).  In  differential  equation  form  this  is 


internal  energy  work 

dt  dt 


enthalpy 


Qi 


(5) 


A differential  equation  for  pressure  can  be  derived  by  adding  the  upper  and  lower 
layer  versions  of  Eq.(5),  noting  that 


dE, 

dt 


d{cymiTi) 

dt 


(6) 


to  obtain 

f = + (7) 

Differential  equations  for  the  layer  volumes,  energy,  density  and  temperature  can 
be  derived  from  definitions  (1)  to  (4)  using  the  quotient  and  product  rules.  These 
differential  equations  are  derived  in  [9]  and  are  summarized  in  Table  1.  Notice  that 
a ^ term  occurs  in  all  but  the  mass  equations.  Handling  the  pressure  equation 
properly  is  then  crucial  for  solving  the  zone  fire  model  both  correctly  and  efficiently. 

In  the  following  sections,  we  examine  the  pressure  equations  which  arise  in  some 
simple  examples  of  single  and  multi-room  zone-fire  models.  The  appropriate  nondi- 
mensionalization  of  the  pressure  equations  helps  to  reveal  the  analytical  character  of 
the  equations  and  allows  us  to  use  analytical  techniques  to  examine  the  nature  of  the 
stiffness  that  arises  from  the  pressure  equations. 
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Table  1:  Conservative  Zone  Modeling  Differential  Equations 


Equation  Type 

Differential  Equation 

mass  of  layer  i 

pressure 

^ + gu) 

energy  of  layer  i 

volume  of  layer  i 

1 

1-H 

1 

II 

density  of  layer  i 

dt  ~ CpT.V,  ((?*  Cj,ThtT,)  ^^  ) 

temperature  of  layer  i 

dt  ~ Cpp.V,  Cp^tT,)  + U ^<  ) 

3 Pressure  Equation  One-Room  Model 

In  Section  3.3.1  of  [9],  the  first  example  is  presented  to  illustrate  the  nature  of  the 
problems  encountered  with  the  integration  of  the  pressure  in  zone  fire  modeling.  This 
example  is  illustrated  in  Figure  1 (Figure  3 of  [9])  and  is  for  a fire  in  a single  room 
with  a small  leak  near  the  floor.  The  major  approximations  made  in  this  model  are 
that  the  heated  gases  never  exit  through  the  vent  (i.e.  the  vent  is  near  the  floor,  while 
the  heated  gas  layer  stays  above  the  leak)  and  that  there  is  no  heat  transfer  to  the 
walls;  i.e.  the  enclosure  hcis  adiabatic  boundaries.  Under  these  conditions,  the  mass 
loss  through  the  vent  equals  a constant  times  the  square  root  of  the  pressure  difference 
across  the  vent,  and  the  enthalpy  loss  through  the  vent  is  a constant  times  the  mass 
loss.  The  starting  equation  for  this  model  is  Eq.  (13)  of  [9]  which  is  equivalent  to 
Eq.  (7),  derived  earlier.  We  repeat  the  analysis  given  there  in  nondimensional  form 
and  with  a slightly  different  notation  because  the  nondimensionalization  gives  some 
additional  insight  and  because  this  example,  which  can  be  solved  analytically,  forms 
the  basis  for  analysis  of  more  complicated  cases. 

dp 

di 

Here  the  notation  is  as  defined  in  [9]  except  that  p is  the  pressure  in  the  room  relative 
to  the  ambient  pressure  outside  the  room  p = p^nc  — Patmi  and  t is  the  time,  qjire  is 
the  constant  fire  (heat)  source,  and  9„ent  is  the  enthalpy  loss  through  the  vent,  defined 


ire  ^vent) 
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Figure  1:  One  Room  Test  Case  Configuration 


as 


^vent  — CpTuent^^ent-^vent  Sgn 


vent 


The  initial  conditions  are  that  p = 0 at  t = 0. 

Now  use  the  following  quantities  with  which  to  nondimensionalize: 


Poo 

T 


(IJi 


CpTygYit  Cvent  -^vent  y/^Pvent 
Vpoo  V 


(IjlT 


(7  - l)9/:re  (7  ^Cp'I^vent^vent-^veni  v>  vent ) 


(8) 


The  quantity  poo  = P final  — Patm  is  the  asymptotic  pressure  rise  in  the  enclosure 
relative  to  the  ambient  pressure  due  to  the  specified  heating  with  the  enclosure  leak 
specified.  Similarly,  the  quantity  r is  the  time  scale  over  which  the  pressure  rise 
occurs.  These  are  the  proper  scaling  parameters  with  which  to  make  the  dependent 
and  independent  variables  nondimensional: 

_P_ 

Poo 

i 

T 


P = 

t = 


6 


p 


Figure  2:  Solution  of  a Non-Dimensionalized  Pressure  Differential  Equation  for  Initial 
Conditions:  p(0)  = —2,  p(0)  = 0 and  p(0)  = 2 

Then  the  equation  for  the  pressure  becomes 


The  solution  to  this  equation  with  the  initial  conditions  that  p = 0 at  f = 0 was  given 
earlier  in  [9]  and  in  [10]. 

The  value  of  the  nondimensionalization  will  be  apparent  in  the  second  example, 
the  two-room  model,  described  below. 

Figure  2 shows  solution  plots  to  Eq.  (9)  for  three  values  of  the  initial  conditions, 
p(0)  = 2.0, 0.0  and  —2.0.  These  plots  were  generated  using  the  software  package 
Mathematica[ll],  particularly  the  command  NDSolve.  The  initial  condition  that 
p = y = 0.0  for  t = X = 0.0  is  the  base  calculation  performed  both  in  [10]  and 
in  [9].  For  all  initial  conditions  shown,  the  solutions  all  converge  at  long  time  to 
p{t  — > cx))  = 1,  the  stable  equilibrium  solution  to  the  problem. 


3.1  Phase  Plane  Analysis  for  One-Room  Model 

The  phase  plane  for  the  one-room  model  is  just  the  line,  but  it  is  instructive  to 
consider  it.  We  will  examine  first  the  case  of  positive  dimensionless  pressure  and 
study  the  solution  for  various  initial  conditions.  Then  in  the  latter  portion  of  this 
section,  we  examine  also  solutions  for  negative  initial  dimensionless  pressures.  In  all 
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cases,  the  solutions  pass,  after  long  time,  to  the  stable  equilibrium  value  of  p = 1, 
which  is  called  the  fixed  point  of  the  equation. 

To  solve  Eq.  (9)  generally  for  positive  initial  pressures,  let  p = and  rewrite 
Eq.  (9) 


Integrating, 


dp 

1 - y/P 


= dt 


2xdx 


I — X 


t — to  = —2 In  \ l — x\  — 2x  C 


where  C is  the  integration  constant,  and  where  x < 0 as  well  as  x > 0. 

Now,  let  the  initial  conditions  be  Penc(fo)  = Penco  at  time  f = fo  = 0,  where 
we  retain  the  symbol  to  for  use  later.  Then  po  = Penc  — Patm  and  po  = Pol  Poo  = 
{Penc  — Patm) /poo-  Now,  since  po  < 0 as  well  as  po  > 0,  and  since  poo  can  be  arbitrarily 
small,  — oo  < Po  < oo.  Hence  the  phase  space  for  the  one-room  model  is  the  whole 
line. 

For  0 < Po  < oo. 


t - to  = + 2{^0  - Vp) 

and  the  phase  diagram  for  either  po  > 1 or  for  po  < 1 is  the  directed  line  segment 
from  Po  to  unity. 

The  solution  to  Eq.  (9)  for  negative  initial  dimensionless  pressures  can  be  found 
as  follows.  For  — oo  < po  < 0, 


Using  p = — x^. 


or,  integrating. 


= dt  = 2 


\ 1 + X 


dx 


f f 0 = 2 In  1 1 + X I — 2x  -f  C' 

where  C is  the  integration  constant,  and  where  x > 0. 


Now,  initially,  at  < = 0,p  = po-  Hence, 


( = 21n-ii^  + 2(\/M-ybi) 
l + VW 

When  p starts  at  negative  values,  it  rises  to  zero  at  t = tc, 

tc  = 2\n 

1 + VN 

For  t > <C5  P > 0,  and  we  use  the  solution  for  p > 0 with  initial  conditions  that  p = 0 
at  t = tc,  namely. 


t — tr  = 2\n 


1 


- 2,/p 


|l-,/p| 

In  this  case,  the  phase  diagram  is  a directed  line  segment  from  po  < 0 to  unity. 


4 Pressure  Equation  Two-Room  Model 


The  second  example,  and  the  more  interesting  one  since  it  illustrates  the  structure  of 
the  mathematical  problem  in  solving  the  pressure  equations,  is  illustrated  in  Figure 
3 (Figure  6 of  [9]).  In  this  example,  there  are  two  rooms.  In  the  first  room,  denoted 
by  subscript  1,  there  is  a fire  and  two  vents,  one  to  the  outside,  which  is  denoted  vent 
1,  and  the  second  vent,  denoted  vent  2,  to  the  second  room.  Again,  the  walls  are 
assumed  to  be  adiabatic.  The  equations  for  this  example  are: 


where 


dpi 

di 

dp2 

di 


7 - 1 

7 - 1 
^2 


{Qfire  Qventl  ^vent2^ 
{Qvent2^ 


Qventl  — C'pTyg'fiiiCygYLtl  -^ventl  (Pl  ) 2p|pi  | 

Qvent2  — C pTvent2^vent2  ■^vent2  P2)y^p\Pl  P2 1 


We  define  the  following  scaling  parameters: 


Pool 


Qfire \ 

, CpTygfitl  Cyg'rHl  -Ayefill  \/2p  J 


2 
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Figure  3:  Two  Room  Test  Case  Configuration 


Poo2 

n 


9/t>e A 

CpTrjent2(^ent2-^vent2\/‘^P  j 

FlPool  Fi 


Qfir 


T2  = 


(7  - l)9/i 

^Poo2 


(7  - 1)  {CpT^entl  Cvent\  ^ventX  y/^P? 

F2  Qfire 


(7  ^)9/i’‘e  (7  1)  ( CpTygjif2^vent2  ■^vent2 


We  use  the  pressure  scale  and  the  time  scale  associated  with  vent  one  as  the  basic 
scaling  parameters  with  which  to  make  the  equations  dimensionless.  Define 


Px  = Px/Poox 

P2  = P2/P00X  (10) 

t = i/ri 


Then,  the  ratios  Poox/Poo2  and  ri/r2  appear  in  the  dimensionless  equations,  which  we 
can  write  as  follows 


^ = 1 - sgn(pi)yH  - y^sgn(pi  -p2)\/\px  -P2\  (11) 

dp2  /tiW  . . A : 
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4.1  Phase  Plane  for  Two-Room  Model 

For  the  phase-plane  analysis,  change  variables  to  pi,Ap  = P\  — P2-  Then  these 
equations  become 

^ = 1 - sgn(pi)yH  - y^^sgn(Ap)yiA^ 


dAp 

dt 


= l-sgn(pi)/H- 


Vl  “I"  V2  /tj  V2 


V2 


r2Vi 


sgn  (Ap)y|Ap| 


These  equations  can  be  rewritten  as  follows;  let 


Pi 

Ap 


X 

y 


a = 


b = 


TjV2 
Vt2Vi 
Vi-f  ^2 
V2 


Then  the  equations  above  become 

^ = 1 - sgn(a:)yjzj- asgn(y)yM 

dy 


= 1 - sgn(a:)yi7|  - afesgn(y)yH 


(12) 


These  equations  are  autonomous,  see  [12]  and  [13]  for  example,  and  therefore  can 
be  reduced  to  a single  first-order  nonlinear  ODE  by  eliminating  time.  Dividing  these 
two  equations,  we  get 

^ _ 1 - sgn  " absgn  {y)y/\y\ 

1 - sgn(x)yH  - asgn(y)yi^ 


or 


^ = \- 
dx 


a{b-  l)sgn(y)yi^ 


1 — sgn  (x  )\/W-  asgn(!/)Y'M 


(13) 


Although  a phase  plane  analysis  formally  starts  from  either  of  these  equations,  when 
we  numerically  integrate  in  the  phase  plane,  we  use  the  parametric  form  of  the  equa- 
tions, Eqs.  (12),  and  then  plot  the  results  in  the  phase  plane. 

If  the  two  room  volumes  are  identical  and  the  conditions  in  the  lower  layer  are  the 
same  for  the  example  illustrated  in  Figure  3 then  the  parameters  a and  b are  given 
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Figure  4:  Solution  Plot  of  Equation  12  for  a = 1,  6 = 2 and  Initial  Conditions 
a;(0)  = y(0)  = 0 

by  a = Avent2Mventi  and  6 = 2.  Large  a then  implies  that  the  vent  connecting  the 
two  rooms  is  large  compared  to  the  vent  connecting  the  first  room  to  the  outside. 

In  general,  the  fixed  point  of  the  system  for  {x{t),y{t))  is  determined  from  the 
equations 


and  is  given  by  xq  = l,yo  = 0.  We  note  an  important  difference,  between  the 
equations  for  {x{t),y{t))  and  the  usual  ones  encountered  in  phase-plane  analysis  [12] 
and  [13];  these  equations  are  not  analytic  around  the  fixed  point.  The  fixed  point  is 
a stable  one  as  determined  by  the  numerical  integrations  described  below. 

Equations  (12)  have  been  integrated  using  the  software  package  Mathematica. 
Once  again,  as  in  the  one-room  example,  we  have  used  the  command  NDSolve  for 
this  first-order  nonlinear  ODE  system.  Figure  4 shows  x{t)  and  y{t)  for  a = 1,  6 = 2 
with  initial  conditions  (I.C.)  a:(0)  = y(0)  = 0.  The  dimensionless  pressure  in  room 
1,  x(f),  starts  at  zero  and  increases  monotonically  to  unity  over  of  order  ten  dimen- 
sionless time  units.  The  pressure  difference  between  room  1 and  room  2,  y{t),  starts 
at  zero,  increases  to  a maximum  of  about  0.1  at  about  one  dimensionless  time  unit 
and  then  decreases  to  zero  again.  The  solutions  displayed  are  well  behaved,  and  the 
numerical  calculation  of  them  encounters  no  particular  difficulty.  However,  when  left 
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x,y 


Figure  5:  Solution  Plot  of  Equation  12  for  a = 4,  6 = 2 and  Initial  Conditions 
x(0)  = y(0)  = 0 


in  dimensional  form  and  coupled  with  the  equations  describing  the  other  properties 
in  the  zone  model,  the  equations  are  stilf.  Therefore,  in  general,  computation  of  solu- 
tions to  these  equations  could  be  difficult  because  of  the  short  pressure-equilibration 
time  scale  compared  to  the  time  scale  for  change  of  the  other  properties,  such  as 
temperature,  layer  height,  etc. 

Figure  5 shows  x{t)  and  y{t)  for  a = 4,6  = 2 with  initial  conditions  (I.C.) 
x(0)  = y(0)  = 0.  The  primary  difference  between  these  plots  and  those  of  Fig. 
4 is  that  the  solution  for  y rises  more  rapidly  as  a function  of  t from  its  initial 
condition  to  a smaller  maximum  of  about  0.01  and  then  decays  to  zero.  Large  a 
is  a condition  on  the  ratio  of  time  scales  and  volumes  of  the  two  rooms.  As  we 
shall  see  in  the  following  subsection,  large  a implies  that  the  equations  are  stiff  and 
a singular  perturbation  analysis  of  the  problem  is  applicable.  We  emphasize  that 
this  problem,  with  a relatively  large  value  of  a,  would  most  probably  cause  difficulty 
using  a numerical  solver  that  did  not  account  for  stiffness  of  the  equations.  As  noted 
above,  when  the  pressure  equations  are  coupled  with  the  equations  for  the  other 
properties  in  a zone  model,  they  become  stiff  due  to  the  disparity  in  time  scales. 
Therefore,  integration  of  these  equations  would  most  likely  cause  even  more  difficulty 
than  integration  of  the  previous  case  illustrated  in  Fig.  4. 

Figure  6 shows  x{t)  and  y{t)  for  a = 0.1,6  = 2 with  initial  conditions  (I.C.) 
a:(0)  = y(0)  = 0.  The  curves  show  a much  more  gentle  time  dependence  than  that 
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Figure  6:  Solution  Plot  of  Equation  12  for  a = 0.1,  6 = 2 and  Initial  Conditions 

x(0)  = t/(0)  = 0 

displayed  in  Fig.  5. 

Figure  7 shows  a phase  plane  plot  of  the  solution  to  Eqs.  (12)  with  a = 1, 6 = 2. 
This  phase  plane  plot  demonstrates  that  a:  = l,y  = 0isa  stable  fixed  point  of  the 
solution  since  all  solutions  progress  toward  this  point  as  time  increases.  The  plot  was 
prepared  by  integrating  Eqs.  (12)  for  thirteen  different  initial  conditions  and  then 
plotting  each  curve  parametrically.  All  solutions  reach  the  stable  fixed  point. 

Figure  8 shows  a phase  plane  plot  of  the  solution  to  Eqs.  (12)  with  a = 4,b  = 2. 
This  figure  shows  that  the  trajectories  of  the  solution  have  become  much  more  angular 
with  nearly  45-degree  lines  joined  to  sections  of  the  x-axis.  This  very  abrupt  behavior 
is  an  indication  that  the  Eqs.  (12)  are  becoming  stiff  for  the  parameters  chosen. 

Figure  9 shows  a phase  plane  plot  of  the  solution  to  Eqs.  (12)  with  a = 0.1,  6 = 2. 
This  figure  shows  that  the  trajectories  of  the  solution  have  become  much  smoother 
than  those  shown  in  either  Figs.  7 or  8. 


4.2  Zero-Order  Singular-Perturbation  Analysis 

Return  to  Eqs.  (11)  and  consider  the  cases  when  the  ratio  of  time  scales  Ti  / T2  becomes 
either  large  or  small;  these  cases  are  actually  the  cases  of  interest  because  physical 
parameters  dictate  that  the  time  scale  ratio  will  often  be  large  or  small,  and  these 
large  or  small  numbers  cause  stiffness  in  the  equations.  We  will  not  perform  a formal 
singular  perturbation  analysis,  but  only  show  how  the  zero-order  behavior  of  the 
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Figure  7:  Phase  Plane  Plot  of  Equation  12fora  = l,6  = 2 
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Figure  8:  Phase  Plane  Plot  of  Equation  12  for  a = 4,  6 = 2 
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Figure  9:  Phase  Plane  Plot  of  Equation  12  for  a = 0.1,  b = 2 


system  can  be  determined  when  the  ratio  of  time  scales  is  large  or  small. 

First,  consider  the  case  when  Ti/r2  ^ 1.  This  will  occur,  for  example,  when 
the  area  of  vent  two,  that  connecting  the  two  rooms,  is  moderately  large,  while  vent 
one,  the  vent  to  the  outside,  is  a small  leak.  When  this  is  the  case,  then,  since  Ti/t2 
multiplies  the  difference  in  room  pressures,  pi  —p2  — »•  0.  If  we  eliminate  (ti  / T2){pi  — P2) 
between  the  two  equations,  we  get 


If  we  now  say  p2  ~ Pi , then 


If  we  choose  the  proper  pressure  and  time  scales,  this  becomes  the  same  as  Eq.  (9) 
for  the  single  room  with  a fire  and  a leak,  but  now  with  a volume  Vi  + W,  the  volume 
of  the  two  rooms. 

Similarly,  when  ri/r2  <C  1,  we  have  the  case  where  the  vent  area  between  the  two 
rooms  is  small  relative  to  the  vent  area  to  the  outside  for  example.  In  this  case,  we 
concentrate  on  the  first  of  Eqs.  (11)  and  note  that  the  term  proportional  to  ri/T2, 
the  term  representing  the  effect  of  the  second  vent,  is  negligible.  Then,  this  equation 
becomes 
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the  equation  for  the  single-room  case  again. 


5 Conclusions 

The  simple  problems  examined  in  this  paper  illustrate  the  nature  of  the  difficul- 
ties long  encountered  when  numerically  integrating  zone  fire  models.  The  pressure 
equations  equilibrate  very  rapidly  compared  to  the  equations  governing  the  other 
dependent  variables  in  the  zone  fire  models.  When  equations  of  this  nature  are  en- 
countered, they  are  referred  to  as  stiff.  The  simple  problems  analyzed  here  illustrate 
the  nature  of  the  stiffness  and  demonstrate  that  proper  nondimensionalization  to- 
gether with  singular  perturbation  analysis  can  provide  insight  into  the  behavior  of 
the  system  for  parameters  of  interest. 

The  methods  can  be  used  to  examine  much  more  general  problems.  For  example, 
two  rooms  connected  with  each  other  and  with  the  outside  in  different  fashion  can 
be  analyzed  similarly  to  the  two-room  example  presented  here.  In  addition,  some 
multiroom  enclosures  have  also  been  analyzed  using  the  nondimensionalization  and 
singular  perturbation  methods  described  herein.  In  the  limit  of  various  leak  sizes 
between  rooms  (or  time  scales  determined  by  the  heat  source,  room  volume  and  leak 
rate),  the  equations  can  be  shown  to  reduce  to  the  one- room  equation  with  redefined 
leak  rates  and  room  volumes,  as  was  done  in  the  two-room  Ccise  illustrated  above.  The 
methods  should  provide  an  opportunity  to  analyze  difficulties  with  stiffness  which  are 
encountered  in  more  general  zone  fire  models. 
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A Zone-Fire  Model  Stiffness 


Stiffness  is  a term  used  to  categorize  a class  of  mathematical  problems  arising  from 
physical  systems  that  possess  multiple  time  scales.  Paradoxically,  solutions  to  stiff 
problems  often  appear  to  change  slowly  yet  have  enormous  computational  require- 
ments when  solved  using  standard  ie  non-stiff  solvers.  These  solvers  are  inefficient 
because  they  do  not  adequately  exploit  the  key  property  of  stiff  systems.  In  stiff 
systems,  the  short  time-scale  phenomena  approach  a quasi-steady  state  rapidly  while 
the  other  phenomena  evolve  on  a much  longer  time  scale. 

The  derivation  of  the  differential  equations  for  a simple  zone  fire  model  are  pre- 
sented here  to  illustrates  the  presense  of  these  multiple  time  scales.  Once  the  stiffness 
is  understood,  analytical  techniques  can  be  used  to  simplify  the  physical  model  into 
a non-stiff  form  or  alternatively,  special  numerical  procedures  can  be  designed  that 
exploit  the  properties  of  stiff  models. 


A.l  A Derivation  Illustrating  the  Presense  of  Stiffness 


We  start  with  the  equations  of  conservation  of  energy  and  mass  for  each  of  the  two 
layers.  These  equations  are  coupled  with  the  state  equations  and  the  various  defini- 
tions. Here,  tildes  are  used  to  denote  dimensional  quantities.  The  notation  is  that 
used  the  main  body  of  the  paper. 

The  energy  equations  are 


dEu 

, -dVu 

— 9flre  "b  9plume  9walls 

(15) 

di 

di 

dh 

, .dVi 

. (16) 

di 

di 

— 9vent  9plume 

The  mass  equations  are 

drhu 

di 

^^fire  “b  ^^plume 

(17) 

II 

i 

TTlvent  ^plume 

(18) 

The  state  equations  are  for  an  ideal  gas  with  constant  specific-heat  coefficients 


p{i)  = Rpu{i)fu{i)  (19) 

p{i)  = RpL{i)TL{i)  (20) 

Eu{i)  = Cvmu{i)fu{i)  = -Vu{i)p{i)  (21) 

7 - ^ 

Eiii)  = CvrhL{i)fi{i)  = -VL{i)p{i)  (22) 

7 - 1 
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The  definitions  are 


II 

mu{i) 

Mi) 

(23) 

p{i)  = 

miii) 

mi) 

(24) 

Vo  = 

= constant 

(25) 

9 vent  — 

CpT]_,(t^77ly^jll 

(26) 

^vent  — 

C^vent-Avenfy/^PLCP  Po) 

(27) 

9plume  — 

CpTj_,(t^7Tlp\yiaie 

(28) 

rUpimne  — 

(29) 

9walls  — 

[S  + 2(L  + W)(H-  ■z)]K(fu  - To) 

(30) 

^fire  is  assumed  to  be  a constant  fire  (heat)  source,  ^vent  is  the  enthalpy  loss  through 
the  vent,  ^piume  is  the  enthalpy  pumped  from  the  lower  layer  into  the  upper  layer  by 
the  plume,  and  ^waiis  is  the  heat  transfer  rate  to  the  walls,  mfireis  the  mass  added  by 
the  fire,  rhventis  the  mass  loss  through  the  vent  and  mpiumeis  the  mass  pumped  by  the 
fire  plume  from  the  lower  layer  into  the  upper  layer.  Q*  is  the  fire  input  parameter 
defined  by  Zukoski  [5] 


9fire 

PoCr,Tos/PlH^ 

and  po,  To  are  reference  density  and  temperature,  g is  the  acceleration  of  gravity  and 
H is  the  height  of  the  enclosure.  S = LW  is  the  floor  area  of  the  enclosure,  where 
L is  the  length  of  the  enclosure  and  W is  its  width.  K is  a heat  transfer  coefficient 
for  heat  transfer  from  the  gas  to  the  enclosure  walls.  Then  Vq  = LWH.  The  control 
volumes  are  around  the  upper  and  the  lower  layers,  and  the  sum  of  the  two  control 
volumes  is  the  total  volume  of  the  enclosure,  which  is  constant. 

An  equation  for  the  pressure  is  found  by  adding  the  energy  equations  for  the  two 
layers,  taking  account  of  the  equation  of  state  for  the  ideal  gas. 


dEu 

dt 


dEi  . d ~ T/\- 

T “b  — 9fire  Qvent  9walls 


T — (T  l)(9ftre  9vent  9walls) 

at 


dp 

di 


1 - 1 
Vo 


(^fire  Q\ent  ^walls) 
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Here  the  notation  is  as  defined  earlier  except  that  p is  the  pressure  in  the  room  relative 
to  the  ambient  pressure  outside  the  room  p = penc  ~ Patm?  S'lid  t is  the  time. 

The  upper-layer  energy  equation  and  the  upper-layer  mass  equation  combine  to 
give  an  equation  for  the  upper-layer  temperature.  The  upper-layer  energy  equation 
becomes: 


1 d , ^ , .dVv 
1 dp  , JVu^  , JVu 
7 d ~ ~ dp 


The  upper-layer  mass  equation  is 


drhu  _ dpuVu 

dt  di 

^[f.{RfvhVu)  - hVuf-{Rtu)] 


9fire  T 9pl  ume  9 walls 
9fire  T 4pl  ume  9 walls 
Qfire  T 9pl  ume  9walls 

— TTlpJujjje  -j-  TTlf^^ 

— CpTu(^TTlp\  ume  T fiifire) 

— (7p7[/(7fipiume  T fi^ftre) 

— C pTu T fi^fire) 


Combining  the  upper-layer  energy  and  mass  equations 


7 


'Z  ~ d ~ ~ dp 

^ RpU h(/  dF  9plume  ^walls  C'p7(/(rhpiujjie  T fil'fire) 


Now  noting  Eq.(28)  and  that  Cp  = -:^Ri 


CppuVu^ifu)  = + 4fire  “ 9walls  “ Cp{fu  - T’L)mpiume  “ CpfuTUfire 


Finally,  the  lower-layer  mass  equation  is 


dmi  dpiVi 

dt  di  — ^plume  ^vent 

Rewrite  the  three  equations  for  the  pressure,  the  upper  layer  gas  temperature  and 
the  lower  layer  mass. 


(^fire  9vent  9wciUs) 

yo 


dp 

di 
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CpPU^U-^iTu)  = + 9fire  — 9walls  — C'p(7c7  — 7L)mpiujne  — C'pT’t/mfire 

at  at 


dhVL 

di 


— fMplume  ^vent 


Use  much  of  the  notation  and  the  scales  of  Zukoski  to  define  dimensionless  quan- 
tities. If  z is  the  dimensional  layer  height,  then 


2 = z/H 
Vl  = VlIVo  = z 
Vu  = Vu/Vo  = l-z 

Define  the  reference  pressure  po,  temperature  To  and  density  po  to  be  related  by 
po  = RpoTo. 

Pl  = hIpQ 
Pu  = PuIpo 
Tl  = Tl/To 
Tu  = fu/fo 


The  time  scale  defined  by  Zukoski  for  nondimensionalization  is 

<Z  = 4^9(8!  H^) 

The  only  tricky  feature  of  the  current  nondimensionalization  is  that  the  pressure  must 
be  defined  carefully.  Let  p = po  + Ap;  i.e.,  the  pressure  is  the  reference  pressure  plus 
the  overpressure  in  the  enclosure.  The  overpressure  must  be  made  dimensionless  with 
respect  to  poo  where 


Poo 


9fire 


C'p  Tl  C vent -d  vent  v^PL  , 


(32) 


Then 


p = p/po  = 1 + Ap/po  = 1 + (poo/po)Ap  = 1 -f  eAp 

where  e = (poo/po)- 

Rewrite  the  mass  and  heat  sources  and  sinks  as  follows: 

= 


TTTfij-g 


(33) 


Qfire 

= poC.fo^H^Q' 

(34) 

^vent 

= poy/gHH^QvyA^/pL^P 

(35) 

9 vent 

— CpToT £,771  ycnt 

(36) 

= CpfopoyJgHH'^Qvy/tTLyJpL^p 

(37) 

^ plume 

(38) 

9plume 

(39) 

QwaHs 

= [S  + 2(L  + W)(H-z)]K(tu-f^) 

(40) 

= - z)]iTu  - 1) 

(41) 

where 

0* 

^ ~ PoC.To^m 

(42) 

^ Cvent-^vent\/2/?oPo 

" PoVgHH^ 

(43) 

1 2 -'d-vent  Csoimd 

(44) 

o 

" CpPoVgHH^ 

(45) 

The  equations 

in  dimensionless  form  then  become: 

dAp 

dt 

l[Q'  - QvVITl\/^P  -Qw[\  + - z)]{Tu  - 

1)1 

II 

i 

7-1^^  dAp,  ^ n , + ^^rr^  ix 

€ + Q (5w[l  + 2 (1  z)\{Tu  1) 

'y  at  b 

- (Q-yl^aZ^I\Tu  - Tl)  - MTu 

dpLVi 

dt 

-{Q'ff^aZ^I^-QvVl^/pL^P 

The  purpose  of  the  derivation  presented  here  is  to  demonstrate  that  the  equations 
describing  a very  simple  case  of  a zone-fire  model  are  stilT.  Zukoski  [5]  has  presented 

numerical  estimates  of  the  magnitudes  of  the  various  parameters  which  appear  in  these 
equations.  In  particular,  he  has  demonstrated  that  leaks  generally  are  large  enough 
in  most  enclosure  fire  scenarios  that  the  overpressure  which  can  develop  is  rather 
small.  (In  fact,  Zukoski  uses  this  fact  to  ignore  any  overpressure  and  make  a quasi- 
steady approximation  for  the  pressure,  assuming  that  it  equilibrates  instantaneously 
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to  the  pressure  outside  the  enclosure,  during  the  enclosure-filling  process.  If  the 
overpressure  were  to  become  significant,  in  most  cases  there  would  be  a structural 
failure  such  as  a window  breaking  to  relieve  this  overpressure.  A possible  exception 
might  be  any  enclosure  which  was  designed  to  accommodate  large  overpressures,  such 
as  a compartment  in  a submarine  for  example.)  When  the  leak  is  large  enough  to 
sustain  only  small  values  of  overpressure,  then  e = ^ <C  1 is  a small  parameter. 
Using  equation  (32)  with  qtu-e  = 100,  OOOIU,  Avent  = Im^,  Cvent  = -68  and  ambient 
density  and  temperature,  e w 10“®  . Since  this  small  parameter  multiplies  the  time 
derivative  of  the  overpressure,  the  system  of  equations  is  stiff,  and  the  culprit  is  the 
overpressure  equation. 

Other  observations  can  be  made  from  this  dimensionless  system  of  equations. 
However,  we  will  note  only  one.  The  state  equation  is  p = 1 + — PlTl  — PuTu- 

When  e = ^ <C  1,  then  p « 1 and  plTl  = ^iPuTjj  = 1-  Then  Pl  ~ 1 and  7T  « 1. 

A. 2 Some  Numerical  Considerations 

As  pointed  out  in  [9]  many  different  (but  analytically  equivalent)  sets  of  differential 
equations  can  be  derived  to  form  a zone  fire  model.  Further,  the  numerical  difficul- 
ties encountered  in  these  models  because  of  stiffness  can  not  be  avoided  simply  by 
exchanging  the  pressure  equation  for  some  other  equation  such  as  temperature,  den- 
sity, or  internal  energy.  As  shown  in  Table  1,  each  zone  modeling  differential  equation 
contains  a ^ term.  For  a one  room  zone  fire  model,  if  the  pressure  is  computed  us- 
ing equation  (8),  the  formula  for  asymptotic  pressure  rise,  and  ^ is  removed  from 
the  other  modeling  differential  equations,  then  the  resulting  approximate  differential 
equations  are  not  stiff  and  a standard  nonstiff  solver  may  be  used.  This  is  essentially 
why  one  room  models  such  as  ASET  [14]  can  use  non-stiff  solvers  in  their  solution. 
For  multi-room  zone  fire  models,  a non-linear  set  of  equations  need  to  be  solved  to 
obtain  the  quasi-steady  state  pressures  for  each  room.  However,  the  class  of  problems 
that  can  be  solved  is  reduced  since  large  pressure  fluctuations  can  not  be  modeled 
properly. 

The  curious  aspect  of  stiff  differential  equations  is  that  the  solution  appears  to  be 
changing  slowly  and  yet  the  computational  costs  of  computing  this  solution  are  enor- 
mous when  using  nonstiff  differential  equation  solvers  such  as  Runge-Kutta  methods. 
The  question  then  is  why  does  it  cost  so  much  to  solve  a problem  whose  solution 
changes  slowly?  To  maintain  stability,  a nonstiff  solver  must  use  a stepsize  that  is 
small  enough  to  track  the  part  of  the  solution  corresponding  to  the  shortest  time 
scale  even  when  this  solution  component  decays  rapidly  to  some  quasi-steady  value. 
This  stepsize  is  much  smaller  than  required  to  accurately  track  the  desired  part  of 
the  solution  which  corresponds  to  one  of  the  longer  time  scales.  So  for  stiff  problems 
the  choice  of  stepsize  is  dominated  by  considerations  of  stability,  not  accuracy. 
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There  is  no  one  definition  of  stiffness  that  is  universally  applied  to  initial  value 
problems.  One  that  is  commonly  applied  is  the  following  (see  [15]).  An  initial  value 
problem,  dyidt  = f{y,t),y{to)  = yo  is  called  stiff  if  the  eigenvalues,  Aj  = Uj  + ivj, 
j = 1, . . . , of  the  Jacobian,  /y,  satisfy 

Wj  < 0,  j = 1, . . . , iV  , 


and 


max  (|u,|)  ^ min  (|u,|)  . 

In  our  case,  the  eigenvalues  with  large  negative  real  parts  correspond  to  the  short 
time  scale  phenomena  (room  pressures). 

B Mathematica  Commands  Used  to  Perform  Cal- 
culations and  to  Produce  Figures 

This  appendix  documents  the  commands  used  to  both  perform  the  calculations  and 
to  produce  the  plots  illustrated  in  this  report.  These  calculations  were  performed 
using  Mathematica[ll].  The  equivalent  analysis  using  traditional  methods  involving 
FORTRAN  would  have  certainly  taken  much  longer;  perhaps  weeks  instead  of  several 
hours.  The  simple  models  for  the  vent  and  fire  used  to  produce  phase  plane  plots  can 
be  made  more  complex  and  hence  more  realistic  by  adding  a few  “equations”  to  the 
following  Mathematica  code. 


***  figure  2 *** 

rhs[y_];  = l.  -Sign[y]  Sqrt[Abs[y]] 

odel  [y0_]  :=NDSolve  [{y  ’ [x]  ==  rhs[y[x]],  y[0]  ==  yO},y  ,{x,0, 10}] 
y2=odel [2] 
yin2=odel  [-2] 
yO=odel [0] 

plot2=Plot  [{Evaluate [y[x]  /.  y2]  .Evaluate [y  [x]  /.  yin2]  .Evaluate [y[x]  /.  yO]}. 
{x  .0 . 10}  .PlotRaiige->{-2.2}.AxesLabel->{"t"  . "p"}] 

♦♦♦♦  Figures  4.5.6  time  dependent  plots  of  x and  y ♦♦♦♦♦ 

ode2 [a_ .b_] := 

NDSolve[{x’ [t]  ==  1 .-Sign[x[t]]  Sqrt  [Abs [x[t]]] -a  Sign[y[t]]  Sqrt  [Abs [y  [t]]]  . 
y’[t]  ==  1 .-Sign[x[t]]  Sqrt  [Abs  [x  [t] ]]  - ab  Sign[y[t]]  Sqrt  [Abs  [y[t]]]  . 
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x[0]  ==  0,y[0]  ==  0},{x,y},-Ct,0,15}] 
ode2[l  ,2] 

plot4=Plot  [{Evaluate  [x[t]  /.  , 

Evaluate[y[t]  /.  '/,]},{t,0,15},AxesLabel->{"t"  ,"x,y">] 
ode2[4,2] 

plot5=Plot  [{Evaluate  [x[t]  /.  */,]  , 

Evaluate[y[t]  /.  '/,]}, {t,0,l},AxesLabel->{"t","x,y"}] 
ode2[.l,2] 

plot6=Plot  [{Evaluate  [x[t]  /.  , 

Evaluate[y [t]  /.  '/,]},{t,0,15},AxesLabel->{"t","x,y"}] 


**♦  Phase  Space  Parametric  Plot  - Fig\ire  7,  a=l,b=2  **♦ 


ode3 [a_ , b_ , x0_ , y0_] : =NDSolve [ 

{x’[t]  ==  1 ,-Sigii[x[t]]  Sqrt  [Abs[x[t]]]-a  Sign[y[t]]  Sqrt  [Abs[y[t]]]  , 
y’[t]  ==  1 . -Sign[x[t]]  Sqrt[Abs[x[t]]]-  ab  Sign[y[t]]  Sqrt  [Abs [y [t]]]  , 
x[0]  ==  x0,y[0]  ==  y0},{x,y>,{t,0,15},MaxSteps->3000] 

xypl=ode3 [1 ,2 , 1 ,2] 
xyp2=ode3 [1 , 2 , 2 , 2] 
xyp3=ode3 [1 , 2 , - . 5 , 1] 
xyp4=ode3 [1 , 2 , 2 , 1] 
xyp5=ode3 [1 ,2, 1 .5,2] 
xyp6=ode3[l,2,-.5, .5] 
xyp7=ode3 [1 ,2 , .5,2] 
xyp8=ode3 [1 ,2,2, .5] 
xynl=ode3[l ,2, 1 ,-l] 
xyn2=ode3[l ,2,- .5,-1] 
xyn3=ode3 [1 , 2 , . 5 , - 1] 
xyn4=ode3[l ,2, .5,-1] 
xyn5=ode3[l ,2,1 .5,-1] 
xyn6=ode3 [1 , 2 , - . 5 , - 1] 

plot7=ParametricPlot [{Evaluate[{x[t] ,y [t]}  /.  xypl] , 

Evaluate  [{x[t]  ,y  [t]  } /.  xyp2]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyp3]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyp4]  , 


26 


Evaluate  [{x  [t]  ,y[t]}  /.  xyp5]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyp6]  , 

Evaluate [*Cx[t]  ,y[t]}  /.  xyp7]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  X3rnl]  > 

Evaluate  [{x  [t]  , y [t]  } / . xyn2]  , 

Evaluate  [-Cx  [t]  ,y[t]>  /.  xyn3]  , 

Evaluate  [{x  [t]  , y [t]  } / . xyn4]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xynS]  , 

Evaluate  [{x  [t]  , y [t]  } / . xyn6]  } , 

{t  ,0 , IS} ,PlotRange->  All , AxesLabel->-C"x"  , " 

***  Phase  Space  Parametric  Plot  - Figure  8 


xyp2=ode3 [4 , 2 , 1 , 2] 
xyp3=ode3 [4 , 2 , 2 , 2] 
xyp5=ode3 [4 , 2 , 1 . 5 , 2] 
xyp6=ode3 [4 , 2 , 0 , 1] 
xynl=ode3 [4 ,2 , 1 , -1] 
xyn3=ode3 [4 , 2 , . 5 , - 1] 
xyn4=ode3 [4 , 2 , 1 . 5 , - 1] 
xyn5=ode3 [4 , 2 , 2 , - 1] 

plot8=ParametricPlot [ 

{Evaluate  [{x  [t]  ,y[t]}  /.  xyp2]  , 
Evaluate[{x[t]  ,y [t]}  /.  xyp3]  , 
Evaluate[{x[t]  ,y [t]}  /.  xyp5]  , 

Evaluate  [{x  [t]  , y [t]  } / . xyp6]  , 

Evaluate  [{x  [t]  , y [t]  } / . xynl]  , 

Evaluate  [{x  [t]  , y [t]  } / . xyn3]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyn4]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyn5]}, 

{t  ,0, 15} ,PlotRcLnge->  All , AxesLabel->{"x"  , " 


***  Phase  Space  Pairametric  Plot  - Figure  9 

xypl=ode3[. 1 ,2,0,0] 
xyp2=ode3 [ . 1 , 2 , 1 , 2] 
xyp3=ode3 [ . 1 , 2 , 2 , 2] 
xyp4=ode3 [ . 1 , 2 , . 5 , 2] 
xyp5=ode3 [ . 1 , 2 , 1 . 5 , 2] 


y">] 


y">] 
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xyp6=ode3 [ . 1 , 2 , 0 , 1] 
xyp7=ode3 [. 1 ,2,2 , 1] 

xynl=ode3 [ . 1 , 2 , 1 , - 1] 
xyn3=ode3 [ . 1 , 2 , . 5 , - 1] 
xyn4=ode3 [.1,2,1. 5,-1] 
xyn5=ode3 [ . 1 , 2 , 2 , - 1] 

plot9=ParainetricPlot  [ 

{Evaluate [-Cx [t]  ,y[t]}  /.  xypl]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyp2]  , 

Evaluate[{x[t]  ,y[t]}  /.  xyp3]  , 

Evaluat  e [{x  [t ] , y [t]  } / . xyp4]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyp5]  , 

Evaluate  [{x  [t]  , y [t]  } / . xyp6]  , 

Evaluat e [{x [t]  ,y[t]>  /.  xyp7]  , 

Evaluate  [{x  [t]  ,y  [t]  } /.  xynl]  , 

Evaluate  [{x  [t]  ,y[t]}  /.  xyn3]  , 

Evaluate [{x [t]  ,y[t]}  /.  xyn4]  , 

Evaluate [{x [t] , y [t] } / . xynS] } , 

{t ,0 ,30> ,PlotRange->  All ,AxesLabel->{"x" , "y"}] 

Phase  Plane  Properties  - Figure  10 
a = 1 
b = 2 

yl [a_ ,x_] :=(1  - Sign[x] Sqrt [Abs [x]] ) “2/a“2 
plotlO=Plot [{yl  [a,x] ,yl[a  b,x] , 

-yl [a,x] , -yl [a  b ,x] } ,{x , -4 ,4} , AxesLabel->{"x" , "y"}] 

save  graphics  to  an  output  file 

stmp=OpenWrite ["f ig2"] 

Display [stmp ,plot2] 

Close  [stmp] 

Ipsfix  -epsf  < fig2  > fig2.ps 

stmp=OpenWrite ["f ig4"] 

Display [stmp ,plot4] 

Close  [stmp] 

Ipsfix  -epsf  < fig4  > fig4.ps 
stmp=OpenWrite ["f ig5"] 

Display [stmp ,plot5] 
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Close [stmp] 

Ipsfix  -epsf  < fig5  > fig5.ps 
stmp=OpenWrite ["f ig6"] 

Display [stmp ,plot6] 

Close [stmp] 

Ipsfix  -epsf  < fig6  > fig6.ps 
stmp=OpenWrite ["f ig7"] 

Display [stmp ,plot7] 

Close  [stmp] 

Ipsfix  -epsf  < fig7  > fig7.ps 
stmp=OpenWrite ["f ig8"] 

Display [stmp , plots] 

Close  [stmp] 

Ipsfix  -epsf  < figS  > figS.ps 
stmp=OpenWrite ["f ig9"] 

Display [stmp , plots] 

Close [stmp] 

Ipsfix  -epsf  < fig9  > fig9.ps 
stmp=OpenWrite ["f iglO"] 

Display [stmp .plot 10] 

Close  [stmp] 

Ipsfix  -epsf  < figlO  > figl0.ps 

print  graphics  *♦* 

PSPrint [plot2] 

PSPrint [plot4] 

PSPrint [plots] 

PSPrint [plots] 

PSPrint [plot7] 

PSPrint [plots] 

PSPrint [plots] 

PSPrint [plot 10] 
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